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ABSTRACT 


The  objective  of  this  study  was  to  model  and  simulate  progressive  failure 
initiating  from  a  notch  tip  in  a  laminated  fibrous  composite  subjected  to  tensile  in¬ 
plane  loading.  The  micro/macro-level  approach  was  used  for  this  study.  The  micro¬ 
level  analysis  used  the  3-D  unit  cell  model  while  macro-level  analysis  used  the  finite 
element  analysis  technique.  A  cross-ply  laminate  with  double-edge  notches  was 
studied  to  investigate  delamination,  fiber  splitting,  and  transverse  matrix  cracking  in 
the  structure.  Numerical  results  were  compared  to  previous  experimental  work. 
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I. 


INTRODUCTION 


A.  BACKGROUND 

Throughout  history  technological  development  has  been  closely  related  to  our 
ability  to  use  and  manipulate  materials  to  suit  our  needs.  Many  modem  composite 
materials  offer  design  engineers  desirable  combinations  of  material  properties 
unavailable  in  more  traditional  materials.  Composites  are  often  used  as  structural 
members  because  of  their  improved  strength-to-weight  ratio,  stiffness-to-weight  ratio, 
toughness,  or  corrosion  resistance.  They  can  also  be  engineered  to  provide  a  particular 
combination  of  properties  by  varying  such  factors  as  constituent  materials,  volume 
fraction  of  constituents,  and  fabrication/manufacturing  processes.  The  primary  drawback 
from  most  composites  is  their  higher  costs.  However,  increases  in  component  costs  are 
frequently  offset  by  decreases  in  systemic  costs  due  to  weight  savings. 

B.  PREVIOUS  WORK 

A  number  of  investigations  have  been  conducted  concerning  damage  and  failure 
in  composite  materials.  In  the  first  of  their  four  paper  series  on  Damage  Mechanics  of 
Composite  Materials,  Kortschot  and  Beaumont  conducted  an  experimental  study  to 
establish  a  qualitative  relationship  between  the  notched  strength  and  terminal  damage 
state  of  double-edge-notched  (DEN)  cross-ply  specimens.  They  also  identified  three 
primary  failure  modes  for  these  specimens:  fiber  splits  in  the  0°  plies,  transverse  ply 
matrix  cracks  in  the  90°  plies,  and  triangular  delamination  zones  at  the  0/90  interfaces. 
These  modes  are  shown  schematically  in  Figure  1.  They  also  found  that  the  angles  at  the 
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apex  of  the  delamination  zones  measured  between  5°  and  10°  for  most  specimens.  [Ref. 

1] 


Figure  1  Failure  Modes  in  a  (90/0)s  Laminate— After  [Ref.  1] 

The  finite  element  method  has  been  used  by  many  investigators  to  model  and 
predict  the  response  of  composite  materials  to  various  loads.  Chen  et  al  [Ref.  2]  and 
Hitchings  et  al  [Ref.  3]  both  approach  the  problem  using  fracture  mechanics  assumptions 
and  strain  energy  release  rate  criteria  for  delamination  propagation.  Davidson  et  al  [Ref. 
4]  simplified  that  approach  further  by  employing  a  special  crack-tip  element  in  the 
expected  damage  region.  Reedy  et  al  [Ref.  5]  modeled  delamination  by  using  an  eight- 
noded  hex  constraint  element  to  connect  shell  elements  of  different  laminae.  If  their 
failure  criterion  was  met,  the  connection  was  considered  broken  and  delamination 
considered  to  have  begun  or  grown. 

Brewer  and  Lagace  [Ref.  6]  approached  delamination  from  the  viewpoint  of  the 
design  engineer  using  known  composite  materials  in  structural  applications.  As  such, 
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they  sought  a  quick  and  efficient  method  for  predicting  delamination.  They  proposed  and 
verified  a  Quadratic  Delamination  Criterion  based  on  the  interlaminar  stresses  calculated 
efficiently  in  Kassapoglou  and  Lagace’s  previous  work  [Ref.  7].  The  primary  advantage 
of  their  methods  is  computational  simplicity;  it  enables  the  designer  to  predict 
delamination  initiation  without  performing  a  complete  finite  element  analysis  of  the 
structure. 

Another  approach  to  failure  prediction  using  the  finite  element  method  involves 
using  a  micromechanical  model  to  separately  examine  fiber  and  matrix  stresses  in 
composite  laminae.  Ardic  et  al  [Ref.  8]  conducted  a  study  examining  composites  of  dis¬ 
similar  laminae  using  this  methodology.  Zhu  et  al  [Ref.  9]  developed  their  Direct 
Micromechanics  Method  (DMM)  and  compared  it  with  phenomenological  results  for 
unidirectional  composites. 

C.  OBJECTIVE 

The  objective  of  this  study  was  to  reproduce  numerically  damage  patterns  in  DEN 
cross-ply  continuous  fiber  composites  that  had  been  observed  experimentally  but  which 
are  still  unexplained  analytically.  To  achieve  this  a  micro-/macromechanical  model  was 
applied  in  conjunction  with  appropriate  failure  criterion.  A  parametric  study  was  then 
conducted  to  investigate  the  effects  of  notch  size  and  specimen  width  on  the  delamination 
zone  size. 
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II.  SOLUTION  METHODOLOGY 


A.  ANALYTIC 

1.  Determination  of  Material  Properties 

In  order  to  effectively  use  composite  materials,  designers  must  be  able  to  reliably 
determine  their  properties  including  their  degradation  as  a  result  of  damage.  Structural 
analysis  of  composite  components  is  conducted  using  one  of  two  approaches:  micro-  or 
macromechanical.  In  macromechanical  analysis,  the  component  is  treated  as  a  pseudo- 
homogeneous  body.  The  effective  material  properties  of  the  pseudo-homogeneous  body 
are  either  predicted  through  “smearing”  or  averaging  the  properties  of  its  constituent 
materials,  or  they  are  determined  experimentally.  The  smearing  calculation  is  generally 
based  on  micromechanical  analysis,  but  after  the  effective  properties  have  been 
determined,  there  is  no  further  need  for  the  properties  of  the  constituent  materials. 

Micromechanical  analysis  uses  the  properties  of  the  constituent  materials  and 
their  interaction  to  predict  the  response  of  the  composite  material.  Three  general 
approaches  to  micromechanical  analysis  are:  the  Rule  of  Mixtures  (ROM),  semi- 
empirical  approaches,  and  methods  of  cells.  Certain  assumptions  are  common  to  each  of 
these  methods;  for  continuous  fiber  composites  they  include:  homogeneous,  linearly 
elastic,  isotropic  matrix  material;  perfect  bonding  between  fiber  and  matrix;  and 
regularly  spaced,  perfectly  aligned,  homogeneous  fibers  [Ref.  10], 

The  Rule  of  Mixtures  is  a  rough  tool  for  predicting  composite  properties  through 
volume  averaging;  it  is  based  on  the  strength  of  materials  approach.  Longitudinal 
properties  result  from  the  assumption  that  fibers  and  surrounding  matrix  material 
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experience  the  same  strain.  Transverse  properties  follow  from  the  assumption  that  fiber 
and  matrix  are  subjected  to  the  same  stress.  Experimental  work  has  shown  that  actual 
longitudinal  properties  of  uniaxial  laminae  are  comparable  to  those  calculated  through 
ROM;  calculated  transverse  properties  have  proven  to  be  less  comparable.  Thus,  ROM  is 
useful  as  a  rule  of  thumb  for  estimating  properties  along  the  fiber  orientation. 

The  semi-empirical  approach  is  characterized  by  the  Halpin-Tsai  Equations  which 
relate  composite  properties,  p,  to  matrix  properties,  pm,  as  follows: 

p  J  ±Ml  0) 

pm  1  -rjVf 


where  the  function  r|  is  of  the  form  shown: 


El 


-l 


7  = 


(2) 


and  4  is  determined  experimentally  and  depends  on  a  variety  of  factors,  including  fiber 
volume  fraction,  fiber  geometry,  loading  conditions,  and  which  property,  p,  is  sought.  Vf 
is  the  volume  fraction  of  the  fiber.  [Ref.  1 1] 

The  Rule  of  Mixtures  and  the  Halpin-Tsai  Equations  provide  means  of 
determining  smeared  properties  of  a  composite  lamina  in  the  linear  elastic  range.  To 
progress  beyond  the  linear  elastic  range  into  damage  analysis  and  degradation  of 
materials  as  a  result  of  damage,  another  approach  is  required.  The  micromechanical 
method  of  cells  was  introduced  by  Aboudi  [Ref.  12];  in  it  he  assumed  that  a  composite  is 
a  periodic  array  of  matrix  and  reinforcement  (fiber)  materials.  The  analysis  then 
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concentrates  on  a  representative  unit  cell  with  displacement  and  traction  continuity 
boundary  conditions  imposed  at  the  interfaces.  The  unit  cell  is  composed  of  four 
subcells,  one  consisting  of  reinforcement  and  three  consisting  of  matrix  material. 
Although  it  appears  to  imply  that  fibers  have  square  cross-sections,  because  interfacial 
conditions  are  imposed  on  average  vice  point- wise  bases,  the  geometry  is  not  constrained. 
Pecknold  [Ref  13]  noted  that  Aboudi’s  model  forms  the  basis  for  a  finite  element  model 
and  conducted  an  investigation  of  a  simplified  unit  cell  model.  Kwon  and  his  colleagues 
have  refined  the  unit  cell  model  to  examine  stresses  in  both  fiber  and  matrix  components 
at  the  micromechanical  level  [Ref.  14-20].  This  model  is  particularly  suitable  for 
analyzing  damage,  for  failure  criteria  can  be  applied  to  both  fiber  and  matrix  materials  at 
the  micromechanical  level.  This  [Ref.  20],  then,  is  the  basis  for  this  current  work. 

2.  Micro-Macro  Interaction 

In  this  work,  a  combination  of  micro-  and  macro-mechanical  approaches  is  used. 
Micro-mechanical  analysis  is  used  to  determine  smeared,  effective  composite  properties. 
Finite  element  analysis  uses  the  smeared  properties  in  the  constitutive  equations  to 
calculate  macro-  or  elemental  stresses  and  strains.  These  stresses  and  strains  can  then  be 
used  in  a  structural  analysis  or,  as  in  this  work,  decomposed  into  micro-  or  subcell 
stresses  and  strains.  Each  subcell  consists  solely  of  fiber  or  matrix  material,  so  failure 
criteria  can  then  be  applied  to  the  microstresses  and  microstrains  to  determine  failure 
initiation,  progression,  or  mode.  If  failure  is  found  to  have  occurred,  the  material 
properties  of  the  constituent  materials  can  then  be  degraded  to  reflect  that  and  through  the 
smearing  process  the  degradation  is  carried  through  to  the  next  iteration  of  the  analysis. 
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Figure  2  illustrates  the  transitions  and  exchange  of  information  between  the  two  types  of 
analysis. 


Micro/Macro  Approach 


STRUCTURE 


Assessment  of  Failure,  Residual  strength  &  stiffness 


Figure  2  Micro/macromechanical  Model 

3.  The  Micromechanical  Cell  Model 

Although  micromechanical  models  for  fibrous  composites  usually  use  only  four 
subcells,  eight  subcells  were  used  in  this  work  because  the  computer  code  used  was 
adapted  from  one  developed  for  particulate  composites  which  generally  use  an  eight 
subcell  model.  The  principles  are  the  same,  but  the  fiber  occupies  two  subcells  along  its 
longitudinal  direction.  The  unit  cell  is  a  three  dimensional  solid,  rectangular, 
parallelepiped  broken  down  into  eight  subcells,  two  containing  fiber  material  and  six 
containing  matrix  material.  Taking  advantage  of  symmetry  only  one  quarter  of  the  total 
cell  is  used;  as  shown  in  Figure  3. 
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Figure  3  The  Micromechanical  Unit  Cell 


The  size  of  each  subcell  is  dependent  on  fiber  volume  fraction.  The  stresses  and 
strains  for  each  unit  cell  or  element  are  the  volume  averaged  subcell  stresses  and  strains 
as  shown  in  Equations  (3)  and  (4) 

=\[v/(?l+ °i)+ -44  )K+ °i +(Tl+ 4)+  (i-V^)VJ + °0]  (3) 

i,j=l,3 


-  1 

SiJ~2l 


v/(4 + ~44  X£l +su+4  +4)+  C1  -44 )  (4 + 4)] 

U=l,3 


where  ay  and  sy  are  the  elemental  stresses  and  strains,  a)j  and  e*  are  subcell  (k=l,8) 
stresses  and  strains,  and  Vf  is  the  fiber  volume  fraction.  Subcells  1  and  2  contain  the  fiber 
and  subcells  3-8  contain  matrix  material.  Stress  continuity  conditions  at  subcell 
interfaces  are  shown  in  Equation  (5). 
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_ 2  3  4 

**n  ”  **n  ?  **n  ”  **ii 

_i  _  3  _2  _ 
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(5) 


Individual  subcells  may  experience  different  strains,  but  the  sums  of  parallel 
longitudinal  strains  must  be  equal.  A  similar  compatibility  condition  is  imposed  for 
transverse  strains.  Mathematically  these  conditions  are  shown  in  Equation  (6): 


Similar  expressions  are  used  for  shear  strain  compatibility. 

The  constitutive  equations  for  each  subcell  follow  the  form  of  the  generalized 
Hooke’s  Law  shown  in  Equation  (7): 

cr“  =  E^,Su  i,  j,  k,  1  =  1 ,3  and  a  =  1 ,8  (7) 
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The  assumption  of  perfect  matrix/fiber  bonding  can  be  relaxed  if  the  interface 
between  fiber  and  matrix  subcells  is  modeled  using  an  equivalent  spring  as  shown  in 
Figure  4.  This  changes  the  strain  compatibility  equations  to  those  of  Equation  (8). 
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Figure  4  Micromechanical  Unit  Cell  with  Spring  Interface 


^  i  +£u — 1 +£i  i  — A\  +£i  i — ^  i  +£i  i 


Simultaneous  solution  of  Equations  (3)  through  (7)  produces  smeared  composite 
material  properties  in  terms  of  fiber  and  matrix  material  properties  and  their  relative 
compositions.  Combining  this  model  with  the  Finite  Element  Method  detailed  in  Section 
B,  the  following  sequence  of  quantities  can  be  calculated:  elemental  displacements  — * 
elemental  strains  — >  subcell  strains  — *  subcell  stresses  —*■  elemental  stresses. 


11 


In  the  study  of  progressive  damage,  therefore,  failure  criteria  applied  to  fiber  or 
matrix  material  at  the  micro-level  will  result  in  corresponding  degradation  of  material 
properties  at  that  level.  This  degradation  will  then  be  carried  through  subsequent 
iterations  to  reveal  the  failure  processes  of  the  composite. 

B.  NUMERICAL 

1.  Finite  Element  Method 

The  Finite  Element  Method  of  numerical  analysis  is  used  to  solve  the  three 
dimensional  elasticity  problem:  to  determine  stress  and  strain  states  of  both  fiber  and 
matrix,  and  to  apply  failure  criteria  for  each  element.  The  finite  element  method  converts 
a  system  of  partial  differential  equations  into  a  larger  system  of  algebraic  equations.  It 
initially  solves  for  nodal  displacement  and  further  processing  calculates  elemental  strains. 
Application  of  the  micromechanical  model  is  used  to  determine  subcell  strains  and 
stresses;  these  then  serve  as  entering  arguments  for  evaluating  failure  criteria  for  each 
element. 

The  derivation  of  the  equations  of  equilibrium  for  a  three  dimensional  solid  in  the 
absence  of  body  forces  can  be  found  in  any  solid  mechanics  textbook  [Ref  21].  Those 
equations  are: 
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dax  dzxy 
ck  dy 
dcry  dr^ 
dy  dx 

^+£k 

dz  dy 


dz 
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+  dz  ° 

+^-o 

ck 


(9) 


At  first  glance  there  appear  to  be  nine  independent  variables  in  these  three 
equations;  however,  application  of  the  moment  equilibrium  condition  reduces  this 
number  to  six: 


T 

=  T 

xy 

yx 

T 

=  T 

xz 

zx 

T 

—  r 

yz 

zy 

(10) 


Below  is  the  unit  solid  element  illustrating  the  sign  convention  used  with  stress 
variables. 


Figure  5  Sign  Conventions 
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The  variables  ‘u’,  ‘v’,  and  ‘w’  are  used  to  represent  displacements  in  the  ‘x’,  ‘y’, 
and  ‘z’  directions,  respectively.  Assuming  small  displacements,  strain  is  then  defined  as: 


_  du  _  _  dv  _  _  dw 

dx’  y  dy’  dy 
du  dv  dw  dv  dw  du 

7xy=lfy  +  lk’  7yz=~^+^  7x:=lk+lk 


(11) 


These  strain-displacement  relationships  can  also  be  expressed  in  matrix  form  as 
shown  in  Equation  (12): 


e. 
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yyz 
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dz 
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dy  dx. 

0  L 

dz  dy 

0 

dz 


d_ 

dx. 


u 
K  v 
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(12) 


The  generalized  Hooke’s  Law  relating  stress  to  strain  can  be  expressed  as: 

{a}  =  [D]{s}  (13) 


where  {a}  is  the  6x1  stress  column  vector,  {s}  is  the  6x1  strain  vector  and  [D]  is  a  6x6 
matrix  of  material  properties.  Thus  Equations  (  9)-(13)  combine  as  fifteen  equations  with 
fifteen  unknowns,  including  several  partial  differential  equations. 

The  derivation  of  the  finite  element  equations  from  the  three  dimensional 
elasticity  equations  above  is  based  on  course  notes  and  commonly  available  current 
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textbooks  [Ref.  22],  The  method  of  weighted  residuals  as  modified  by  Galerkin  is 
employed  to  develop  the  displacement  based  finite  element  equations.  Eight  noded 
isoparametric  rectangular  parallelepiped  elements  are  used  for  the  formulation.  Linear 
shape  functions  are  used  for  analytical  and  computational  simplicity.  Detailed  derivation 
according  to  these  precepts  is  continued  below. 

The  first  step  in  converting  the  partial  differential  equations  into  algebraic 
equations  is  to  apply  the  method  of  weighted  residuals  to  the  three  dimensional  stress 
equilibrium  Equations  (9).  To  this  end,  each  of  the  three  equations  of  Equation  (9)  are 
multiplied  by  a  weight  function,  Wj  (i=l,  2,  3),  which  is  continuous  over  the  physical 
domain  of  the  problem.  Then,  the  three  new  product  equations  are  integrated  over  the 
entire  problem  domain.  The  goal  is  to  chose  weight  functions  Wj  which  are  orthogonal  to 
the  initial  residuals  of  the  equilibrium  equations  such  that  the  integral  of  their  product  is 
zero.  If  ‘V’  is  the  domain  volume  of  the  problem,  the  weighted  residual  equations  are 
shown  in  Equation  (14): 


dz  'l 

+  atndv=  0 

+^]jW2dv=0 

+&)w3dV=0 


(14) 


At  this  point  it  should  be  noted  that  boundary  conditions  must  be  specified  for  Equation 
(9).  The  boundary  conditions  may  be  specified  as  either  (a)  essential  or  geometric 
boundary  conditions  in  which  some  surface  displacements  are  specified,  (b)  natural  or 
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stress  boundary  conditions  in  which  surface  tractions  are  specified,  or  (c)  a  combination 
of  these  types  of  boundary  conditions.  Before  applying  boundary  conditions,  further 
manipulation  of  the  weighted  residuals  is  required  as  shown  by  Equation  (15). 
Integration  by  parts  will  be  performed  in  order  to  weaken  the  continuity  requirements  on 
the  approximate  solution. 
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In  Equation  (15)  ‘S’  is  the  domain  surface  boundary,  ‘V’  is  the  domain  volume,  and  ‘nx’, 
‘ny’,  and  ‘nz’  are  outward  unit  normal  directional  cosines  in  the  ‘x\  ‘y\  and  ‘z’ 
directions,  respectively.  The  boundary  stress  conditions  are  defined  by  surface  tractions 
as  shown  in  Equation  (16): 

<t>x  =  <rxnx  +  rxyny  +  r x:nz 

(h=TnJt-(7n+Tn 

Yy  Lxy,lx  w ylly  ^  V  r 

<f>z  =  rx:nx  +  xrny  +  ar.n. 

Equations  (15)  and  (16)  can  be  combined  and  expressed  in  matrix  form  as: 
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Equation  (17)  can  be  further  modified  by  separating  the  column  vector  inside  the  volume 
integral  into  the  product  of  a  3x6  matrix  and  a  6x1  vector.  This  step,  shown  in 
Equation(18),  isolates  the  weight  function  derivatives  in  the  matrix  from  the  stress  terms 
in  the  vector. 
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(18) 


Substitution  of  the  strain-displacement  relation  of  Equation  (11)  into  the 
constitutive  relations  given  by  Equation  (13)  provides  another  useful  identity  as  shown  in 
Equation  (19). 
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Equation  (12)  is  now  substituted  for  the  column  vector  on  the  left  hand  side  of  Equation 
(19).  Equation  (18)  and  the  modified  Equation  (19)  are  substituted  into  Equation  (17)  to 
produce  Equation  (20). 
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(20) 


The  preceding  manipulations  of  the  elasticity  relations  have  isolated  displacement 
as  the  desired  quantity.  Now,  to  discretize  the  problem  for  an  algebraic  computational 
solution,  assume  that  over  a  small  domain  each  of  the  displacements  can  be  represented 
by  a  polynomial.  Discretization  divides  a  three  dimensional  domain  into  small-but-finite 
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volume  elements.  In  this  formulation  each  finite  element  has  eight  nodes;  each  node  can 
have  three  orthogonal  displacements.  Each  element,  therefore  has  a  total  of  twenty-four 
degrees  of  freedom  (dof).  Calling  ‘uj’,  ‘v*’,  and  ‘w;’  the  displacements  at  each  node, 
elemental  displacements  can  be  defined  in  terms  of  polynomial  shape  functions  as 
follows: 

8  8  8  ni'i 

u=YJHiui;  v  =  £#>,.;  w  =  YJHiwi 


The  first  partial  derivatives  of  the  shape  functions  exist  as  shown  in  Equations  (22). 
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The  3x1  elemental  displacement  vector  can  be  expressed  as  the  product  of  a  3x24  shape 
function  matrix  and  a  24x1  nodal  displacement  vector.  The  product  of  the  6x3  partial 
differential  matrix  of  Equation  (20)  and  the  3x24  shape  function  matrix  are  typically 
combined  into  one  6x24  matrix.  This  matrix  is  referred  to  as  the  ‘B’  matrix  in  this 
formulation.  In  the  shorthand  notation  shown  in  Equation  (23),  the  ‘B’  matrix  can  be 
partitioned  into  sub-matrices,  ‘Bi’,  where  i=l  to  8. 

W=p,][^][^][^][4]MW[^]l  (23> 


In  this  notation,  the  sub-matrices  are  defined  in  Equation  (24). 
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The  Galerkin  method  only  requires  that  the  weighting  function  be  continuous  over  small 
discrete  intervals  which  correspond  to  the  sides  of  the  finite  elements.  Based  on  the 
Galerkin  weighting  functions  shown  above,  the  weighting  function  matrix  of  Equation 
(20)  can  now  be  written  in  terms  of  the  shape  functions: 
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Equation  (26)  can  be  incorporated  in  Equation  (20)  as  follows: 


[dS 


(27) 


In  Equation  (27),  {d}  is  the  displacement  vector  which  has  grown  from  the  3x1 
vector  of ‘u’,  V,  and  ‘w’  to  a  24x1  vector  of  the  nodal  displacements  ‘uj’,  ‘v*’,  and  ‘wj\ 
The  3x1  vector  in  the  surface  integral  on  the  right-hand  side  has  an  8x1  vector  in  each 
term  so  that  the  entire  integrand  is  shorthand  for  a  24x1  vector  of  discretized  surface 
boundary  tractions.  The  original  three  equilibrium  elasticity  partial  differential  equations 
have  now  been  transformed  to  a  matrix  equation  in  24  terms  for  each  solid  element.  The 
integrals  in  Equation  (27)  may  be  easily  solved  if  simple  shape  functions  are  chosen  and 
if  the  modeled  geometry  is  relatively  simple. 

If  the  modeled  geometry  is  complex,  the  finite  element  geometry  can  be 
simplified  by  use  of  a  transformation  mapping  to  another  reference  space.  In  order  to 
map  ‘x’,  ‘y’,  and  ‘z’  coordinates  of  an  irregularly  shaped  element  onto  ‘r’,  ‘s’,  and  ‘t’ 
coordinates  for  a  rectangular  parallelepiped,  a  Jacobian  transformation  matrix  is  required. 
The  Jacobian  transform  matrix  scales  each  component  of  the  original  space  to  a  new 
space.  Equation  (28)  shows  the  Jacobian  for  an  ‘xyz’  system  mapped  to  an  ‘rst’  system: 
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In  terms  of  the  shape  functions  and  nodal  points  the  Jacobian  is  shown  in  Equation  (29). 
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Since  the  finite  element  integral  equation  includes  the  partial  derivatives  of  the  shape 
functions  with  respect  to  the  ‘xyz’  coordinate  system,  the  inverse  of  the  Jacobian  is 
required.  Let  [J]*1=[r],  where  [T]  is  a  3x3  matrix.  The  shape  function  derivatives  with 
respect  to  the  ‘xyz’  system  are  now  expressed  in  Equation  (30)  with  respect  to  the  ‘rst’ 
system. 
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Equation  (30)  is  used  to  calculate  the  [B]  and  [B]T  matrices  in  the  ‘rst’  system.  Equation 
(29)  is  used  to  transform  the  volume  differential.  Application  of  these  transformations  to 
the  left-hand  side  of  Equation  (27)  is  shown  in  Equation  (31). 
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(31) 
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The  transformation  to  the  ‘rst’  coordinate  system  results  in  simplified  finite 
element  integrals.  The  resultant  transformed  elements  are  termed  isoparametric  elements. 
The  term  isoparametric  refers  to  the  equality  between  the  degree  of  the  equations  for 
transforming  ‘xyz’  into  ‘rst’  coordinates  and  the  degree  of  the  shape  functions  for 
estimating  displacement.  The  shape  functions  in  the  ‘rst’  system  are  shown  in  Equation 
(32). 
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Gauss  Quadrature  is  used  to  numerically  calculate  the  integrals.  The  volume 
integral  is  redefined  as  the  triple  summation  from  1  to  the  number  of  integration  points 
(NIP)  of  the  integrand  evaluated  at  Gauss  integration  points  (r„  Sj,  and  tj)  multiplied  by 
weight  factors  as  shown  in  Equation  (33). 
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The  results  of  the  numerical  integration  may  vary  over  the  elements  in  the  domain  of  the 
model.  Each  of  these  results  can  be  expressed  as  a  24x24  elemental  stiffness  matrix, 
[Ke]. 

Recall  the  surface  traction  boundary  conditions  from  the  right-hand  side  of  Equation  (20). 
The  integration  of  the  directional  component  of  the  applied  stress  over  the  element 
surface  area  to  which  the  stress  is  applied  is  equal  to  the  external  force  applied  to  the 
element.  The  result  of  the  integration  is  a  24x1  vector  {Fe}  which  is  equivalent  to  the 
force  in  the  ‘x’,  ‘y\  and  ‘z’  direction  applied  to  a  surface  of  the  solid.  Substituting  these 
resultant  terms  into  Equation  (20)  completes  the  finite  element  derivation  shown  in 
Equation  (34). 

<34> 

The  finite  element  method  involves  combining  these  elemental  matrix  equations 
with  given  boundary  conditions  to  form  a  large  system  of  simultaneous  equations  to  be 
solved  numerically. 

2.  Implementation 

The  actual  finite  element  analysis  for  this  study  was  conducted  using  a 
FORTRAN  driver  for  three-dimensional  static  analysis  based  in  large  part  on  the 
subroutines  developed  by  Akin  [Ref.  23].  Pre-  and  post-processing  were  conducted  using 
various  MATLAB  programs. 
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The  preprocessor  generated  an  appropriately  formatted  text  file  containing 
material  properties,  nodal  coordinates,  element  connectivity,  boundary  conditions,  and 
loading  data. 

The  FORTRAN  program  used  several  subroutines  to  read  the  input  data  and  store 
it  semi-dynamically.  All  data  was  stored  in  one  of  two  massive  arrays,  one  for  integers 
and  one  for  floating  point  variables,  and  pseudo-pointers  were  used  to  track  the  location 
of  the  first  element  of  each  matrix  of  interest.  The  memory  allocation  for  these  two 
arrays  was  calculated  based  upon  the  input  data.  After  the  input  and  housekeeping 
subroutines  were  completed,  the  heart  of  the  analysis  was  controlled  by  subroutine 
FEMST.  It  called  subroutine  EL3D8A  to  calculate  elemental  stiffness  matrices  and 
assemble  them  into  the  system  stiffness  matrix.  The  stiffness  matrices  were  calculated 
using  the  micromechanical  model  and  Gauss  Quadrature  with  two  integration  points  in 
each  coordinate  direction  for  a  total  of  eight  integration  points  per  element.  Subroutine 
MA3D6F  calculated  material  property  matrices  for  composite  elements,  including 
transforming  those  matrices  to  account  for  fiber  orientation.  After  the  system  matrix  was 
assembled,  boundary  conditions  were  applied,  and  Equation  (34)  was  solved  for  the  nodal 
displacements  using  subroutines  FACTOR  and  SOLVE.  FACTOR  uses  the  Cholesky 
method  to  factor  the  square  system  stiffness  matrix  into  the  product  of  a  lower  triangular 
matrix  and  its  transpose.  SOLVE  uses  Cholesky-Gauss  methods  of  forward  and  back 
substitutions  to  complete  the  solution. 

In  the  process  of  calculating  the  residual,  elemental  stresses  and  strains  are 
calculated  and  then  decomposed  into  subcell  stresses  and  strains.  These  values  can  then 


25 


be  used  to  evaluate  matrix  and  fiber  failure  criterion  and  degrade  the  appropriate 
corresponding  material  properties  if  necessary. 

The  FORTRAN  program  writes  its  output  to  five  main  files:  OUTPUT,  DISPL, 
MATSTS,  FIBSTS,  and  STRESS.  OUTPUT  is  essentially  an  echo  of  the  input  file; 
DISPL  is  the  u,  v,  and  w  displacements  of  each  node;  MATSTS  contains  the  matrix 
subcell  stresses  at  each  integration  point;  FIBSTS  contains  the  fiber  subcell  stresses  at 
each  integration  point;  and  STRESS  contains  the  complete  stress  state  of  each  element. 
These  data  files  were  then  read  into  various  MATLAB  programs  for  further  analysis  and 
visualization. 

C.  FAILURE  CRITERION 

The  nature  of  the  failure  modes  being  investigated  in  this  study,  delamination  and 
splitting,  dictates  that  matrix  material  behavior  be  of  primary  interest.  These  phenomena 
are  of  interest  because  they  can  cause  a  structure  to  fail  at  much  lower  than  expected 
loads.  Delamination  is  an  inter-laminar  process;  as  such,  the  inter-laminar,  through¬ 
thickness,  or  z-direction  stresses  are  expected  to  be  the  controlling  factor  in  its  initiation. 
Brewer  and  Lagace  [Ref.  6]  proposed  the  following  criterion  for  delamination: 
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(35) 


Where  Zs  is  the  shear  strength,  Z*  and  Zc  are  normal  strengths  in  tension  and 

« 

compression,  respectively;  and  xxz  and  xyz  are  interlaminar  shear  stresses,  a1  is 
interlaminar  tensile  stress,  and  crc  is  interlaminar  compressive  stress.  This  criterion  will 
be  the  starting  point  of  this  investigation  in  which  the  failure  criterion  will  be  applied  to 
generate  a  delamination  zone  similar  to  that  observed  experimentally.  Assuming  that  the 
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isotropic  matrix  material  is  twice  as  strong  in  tension  as  it  is  in  shear  and  four  times  as 
strong  in  compression  as  it  is  in  tension,  Equation  (35)  can  be  rearranged  in  terms  of  a 
single  strength  value  as  shown  in  Equation  (36): 


(36) 


In  this  form,  this  criterion  can  be  evaluated  without  necessarily  knowing  the 
strength  of  the  material  in  question.  A  stress  profile  calculated  from  the  left-hand  side  of 
Equation  (36)  can  be  plotted  and  used  to  evaluate  the  suitability  of  the  criterion. 

Similar  criterion  was  applied  using  stresses  transverse  to  the  fiber  direction  to 
evaluate  fiber  splitting  and  transverse  cracking. 
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III.  RESULTS  AND  DISCUSSION 


A.  COMPARISON  OF  2-LAYER  MODEL  TO  3-LAYER  MODEL  TO 

EXPERIMENTAL  RESULTS 

In  their  experimental  work  Kortschot  and  Beaumont  [Ref.  1]  applied  tensile 
loading  to  double-edge-notched  graphite-epoxy  specimens  of  various  cross-ply  lay-ups 
and  sizes.  They  observed  that  delamination  along  the  0/90  interfaces  progressed  from  the 
notch  tip  in  a  triangular  zone  upward  and  toward  the  center  of  the  specimen.  As  the  splits 
in  0°  plies  grew,  the  aspect  ratio  of  the  delamination  triangle  remained  constant.  The 
measured  angles  at  the  apex  of  each  specimen’s  delamination  triangle  were  found  to 
consistently  remain  between  5°  and  10°. 

In  this  study  two  models  were  used  to  simulate  the  behavior  observed  by 
Kortschot  and  Beaumont.  A  two-layer  model  consisting  of  two  laminae,  one  with  fibers 
oriented  along  the  direction  of  loading  (0°)  and  one  with  fibers  oriented  transverse  to  the 
loading  direction  (90°),  was  compared  with  a  three-layer  model  consisting  of  the  same 
two  laminae  plus  a  thin  “delamination  layer”  consisting  solely  of  isotropic  matrix 
material  inserted  between  the  two  laminae.  The  same  material  properties  and  geometry 
were  used  for  each  model.  The  material  properties  of  the  two  laminae  of  each  model 
differed  only  by  their  orientation.  The  properties  of  the  constituents  of  the  specific 
composite  used  in  the  experimental  work  were  unavailable,  but  values  were  chosen  so 
that  they  were  within  approximate  ranges  for  graphite  fibers  and  epoxy  matrix  and  that 
their  smeared  values  were  comparable  to  those  cited  by  those  investigators;  the  values 
used  are  shown  in  Table  1.  The  lay-up  modeled  was  (90/0)s,  that  is,  two  0°  plies 
sandwiched  between  a  pair  of  90°  plies.  Geometric  and  laminar  symmetry  conditions 
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allowed  for  a  proper  model  of  the  specimens  using  only  one-eighth  of  the  original. 
Planes  of  symmetry  were  the  mid-plane  both  longitudinally  and  transversely  and  the 
plane  between  the  0°  plies.  Nodes  located  on  these  planes  were  constrained  in  the 
direction(s)  normal  to  their  respective  planes  of  symmetry.  As  in  the  experimental  work, 
the  modeled  notch  angle  was  60°. 


Table  1  Lamina  Dimensions  and  Material  Properties 


Width 

5mm 

Length 

10mm 

Thickness 

1.5mm 

Notch  Length/Specimen  Width  (2a/w) 

0.125 

Fiber  Longitudinal  Young’s  Modulus  (Efi) 

202.5  GPa 

Matrix  Young’s  Modulus  (Em) 

3.45  GPa 

Fiber  Transverse  Young’s  Modulus  (Eft) 

15.75  GPa 

Fiber  Poisson  Ratio  (12  direction) 

0.29 

Fiber  Poisson  Ratio  (23  direction) 

0.25 

Matrix  Poisson  Ratio 

0.35 

Fiber  Shear  Modulus  (12  direction) 

15.75  GPa 

Volume  Fraction 

0.66 

To  evaluate  the  suitability  of  the  two-layer  model,  the  left-hand  side  of  Equation 
(36)  was  applied  to  the  volume  averaged  matrix  subcell  stresses  for  each  element  in  the 
0°  lamina.  This  profile  is  shown  for  the  entire  modeled  section  in  Figure  6  and  for  the 
refined  mesh  area  in  Figure  7. 
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Figure  6  Delamination  Zone  0°  Ply  2  Layer  Model 
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Figure  7  Delamination  Zone  0°  Ply  2  Layer  Model 
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The  same  sort  of  profile  was  generated  for  the  three  layer  model  based  on  the 
stress  states  of  the  elements  in  the  delamination  layer.  These  are  shown  in  Figure  8  and 
Figure  9. 


Figure  8  Delamination  Zone  3  Layer  Model 
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Figure  9  Delamination  Zone  3  Layer  Model 


It  is  apparent  from  Figures  6-9  that  the  three-layer  model  more  closely  resembles 
the  delamination  zone  observed  experimentally.  The  apex  angle  of  the  delamination 
triangle  (superimposed  in  Figure  9)  in  the  three  layer  model  is  approximately  10°, 
comparing  favorably  with  experimental  observations.  While  the  two  layer  model’s  stress 
profile  does  not  have  the  characteristic  triangular  shape  of  a  delamination  region,  a 
similar  profile  generated  from  its  transverse  matrix  subcell  stresses  does  compare 
favorably  with  observed  0°  ply  splitting.  This  profile  is  shown  in  Figures  10  and  1 1 . 
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Figure  10  Fiber  Splitting  in  0°  Ply,  2  Layer  Model 
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Figure  1 1  Fiber  Splitting  in  0°  Ply,  2  Layer  Model 


The  three  layer  model  not  only  displays  a  more  realistic  delamination  zone,  the 
transverse  stress  profiles  in  the  0°  and  90°  laminae  also  reflect  fiber  splitting  and 
transverse  cracking,  respectively.  As  observed  experimentally  and  expected  with 
graphite  fibers,  transverse  cracking  permeates  the  90°  lamina.  These  effects  are 
illustrated  in  Figures  12-14. 
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Figure  12  Fiber  Splitting  in  0°  Ply,  3  Layer  Model 
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Figure  13  Fiber  Splitting  in  0°  Ply,  3  Layer  Model 
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B.  PARAMETRIC  STUDY 


A  parametric  study  varying  notch  length  and  specimen  width  was  conducted  to 
further  evaluate  the  quadratic  delamination  criterion  as  applied  to  the  three  layer  model. 
Using  the  original  crack  length  to  width  ratio  (2a/w=0.125),  samples  5,  10,  20,  and  40 
mm  wide  were  simulated.  The  resulting  delamination  zones  are  shown  in  Figures  15-18. 
The  apex  angles  of  those  zones  ranged  from  6°-12°. 

In  the  pre-processor  used,  the  nodes  were  placed  as  a  function  of  crack  length  and 
specimen  width  and  length;  as  such,  each  of  these  different  specimens  has  the  same 
number  of  elements.  Therefore,  the  wider  specimens  have  larger  elements  and  some 
detail  is  lost.  Nevertheless,  the  influence  of  the  fiber  splits  and  the  general  triangular 
shape  of  the  delamination  zones  can  be  seen. 
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Figure  15  Delamination  Zone  w=5mm 
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Figure  16  Delamination  Zone  w=10mm 
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Figure  17  Delamination  Zone  w=20mm 
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Delamination  Zone  2a/w=0.125 
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Figure  18  Delamination  Zone  w=40mm 


For  the  same  four  widths  crack  length  to  width  ratio  was  also  varied.  Figures  19- 
22  show  delamination  zones  for  those  specimens.  For  the  first  three  of  these  specimens 
apex  angles  range  from  6°-13°.  Once  again  the  mesh  effects  appear  to  suggest  that 
delamination  is  a  size-dependent  phenomenon.  This  coarsening  is  most  apparent  in  the 
40mm  specimen,  Figure  22. 
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Figure  19  Delamination  Zone  w=5mm 


Delamination  Zone  2a/w=0.25 


Figure  20  Delamination  Zone  w=  10mm 
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Delamination  Zone  2a/w=0.25 


IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  combined  micro-  and  macromechanical  approach  to  analysis  of  composite 
materials  is  an  effective  method  of  exploring  and  simulating  the  initiation  of  damage.  It 
efficiently  enables  the  incorporation  of  constituent  homogeneous  material  behavior  into 
analysis  of  heterogeneous  materials.  In  this  work  micro-/macromechanical  analysis  was 
used  to  simulate  each  of  the  three  experimentally  observed  modes  of  failure  for  DEN 
cross-ply  fibrous  composites.  A  thin  layer  of  homogeneous  matrix  elements  was 
necessary  to  properly  model  observed  delamination  behavior.  The  delamination  layer 
made  the  transmission  of  stresses  and  damage  across  the  laminar  boundaries  conveniently 
observable.  The  quadratic  delamination  criterion  applied  to  delamination  layer  elements 
did  correctly  predict  delamination  zone  shape  for  a  variety  of  geometries.  Proper  mesh 
refinement  is  necessary  to  correctly  simulate  specimen  behavior  in  the  vicinity  of  the 
notch  tip. 

Future  research  in  this  area  should  include  both  experimental  and  simulated 
examinations  of  multi-axial  loading.  This  work  examined  failure  modes  that  occurred 
exclusively  in  the  matrix  regions;  incorporation  of  fiber  failure  and  degradation  should  be 
included  in  future  work. 
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